__global__ void SilhouetteEdgeDetect(unsigned char* surfNorm, unsigned char* surfPos, size_t pitch, uint* AAFlag, uint* AANum, ABSSURF_SHARE_INFO asi, int* PatchIndex, int nPatch, uint pixelBits, uint patchBits)
{
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=cRenderParam.renderWidth || y>=cRenderParam.renderHeight)			return;
	//uint pixelI=x+y*cRenderParam.renderWidth;

	float* pixel = (float*)(surfPos + y*pitch) + 4*x;
	float depth = pixel[3], left, right, up, down;

	left = right = up = down = depth;
	if(x>0)								left  = *(pixel-1);
	if(x<cRenderParam.renderWidth-1)	right = *(pixel+7);
	if(y>0)								up    = *(pixel-pitch/4+3);
	if(y<cRenderParam.renderHeight-1)	down  = *(pixel+pitch/4+3);

	int aaDetected = 0;
	float threshold = 0.15;
	if( (depth==0&&left+right+up+down!=0) || ( left*right*up*down!=0 && ((2*depth-left-right)>=threshold || (2*depth-up-down)>=threshold ) ) )
	{
		aaDetected=1;
		float* pixelN = (float*)(surfNorm + y*pitch) + 4*x + 3;

		float maxD=fabs(depth-left),d1=fabs(depth-right),d2=fabs(depth-up),d3=fabs(depth-down);
		int maxI=0;
		if(d1>maxD)	{maxD=d1;maxI=1;}
		if(d2>maxD) {maxD=d2;maxI=2;}
		if(d3>maxD) {maxD=d3;maxI=3;}
		*pixelN = 1.1+maxI;
		if(maxI==0)	pixel-=4;
		if(maxI==1)	pixel+=4;
		if(maxI==2)	pixel-=pitch/4;
		if(maxI==3)	pixel+=pitch/4;

		//if(fabs(depth-left)>threshold)		{*pixelN = 1.1; pixel-=4;}
		//else if(fabs(depth-right)>threshold)	{*pixelN = 2.1; pixel+=4;}
		//else if(fabs(depth-up)>threshold)		{*pixelN = 3.1; pixel-=pitch/4;}
		//else if(fabs(depth-down)>threshold)	{*pixelN = 4.1; pixel+=pitch/4;}
	}
	float Pos[3] = {pixel[0], pixel[1], pixel[2]};
	

	//if(aaDetected==1)	{pixel[0]=1.0; pixel[1]=0.0; pixel[2]=0.0;}return;

	
	uint thid = threadIdx.x+threadIdx.y*blockDim.x;
	__shared__ uint scanI[SCAN_THREADN];
	scanI[thid] = aaDetected;

	//PreScan***********************************************************************************
	int n = SCAN_THREADN;

	uint* g_idata=scanI, *g_odata=scanI;
	__shared__ uint sumN;
	__shared__ uint preSum;

	__shared__ uint temp[SCAN_THREADN];
	int offset = 1;
	int bankOffsetA = CONFLICT_FREE_OFFSET(thid);
	temp[thid + bankOffsetA] = g_idata[thid];
	for (int d = n>>1; d > 0; d >>= 1) // build sum in place up the tree
	{
		__syncthreads();
		if (thid < d)
		{
			int ai = offset*(2*thid+1)-1;
			int bi = offset*(2*thid+2)-1;
			ai += CONFLICT_FREE_OFFSET(ai);
			bi += CONFLICT_FREE_OFFSET(bi);
			temp[bi] += temp[ai];
		}
		offset *= 2;
	}
	if (thid==0) 
	{ 
		int lastI = n-1 + CONFLICT_FREE_OFFSET(n-1);
		sumN = temp[lastI];
		preSum = atomicAdd(AANum, sumN);
		temp[lastI] = 0; 
	} // clear the last element
	__syncthreads();
	if(sumN==0) return;

	for (int d = 1; d < n; d *= 2) // traverse down tree & build scan
	{
		offset >>= 1;
		__syncthreads();
		if (thid < d)
		{
			int ai = offset*(2*thid+1)-1;
			int bi = offset*(2*thid+2)-1;
			ai += CONFLICT_FREE_OFFSET(ai);
			bi += CONFLICT_FREE_OFFSET(bi);
			float t = temp[ai];
			temp[ai] = temp[bi];
			temp[bi] += t;
		}
	}
	__syncthreads();
	g_odata[thid] = temp[thid + bankOffsetA];
	//PreScan***********************************************************************************


	if(aaDetected)
	{
		uint outInd = preSum+g_odata[thid];		

		uint patInd[3];

		GetBoxInd(patInd, Pos, asi.box0, asi.boxSize);

		uint patchI = patInd[0]*asi.segShift[0]+patInd[1]*asi.segShift[1]+patInd[2]*asi.segShift[2];
		if(PatchIndex[patchI]<0)
		{
			//AAFlag[0] = ((((x<<pixelBits)|y)<<patchBits)|3333);
		}
		else
		{
			patchI = (uint)(PatchIndex[patchI]);
			AAFlag[outInd] = ((((x<<pixelBits)|y)<<patchBits)|patchI);
		}
		//AAFlag[outInd] = (x<<10)|y;
	}

}

__global__ void RSAARayCasting(unsigned char* surface, unsigned char* surfNorm, size_t pitch, int degSum, uint* IntervalPad,	float* RayDir, uint sortPadNum, ABSSURF_SHARE_INFO asi, uint  pixelBits, uint patchBits, int* PatchIndexInverse, float* FuncData)
{
	int thid = threadIdx.x;
	int thidA = thid+__mul24(blockIdx.x, blockDim.x);

	__shared__ uint patchInfo[PIXELTHREADS];
	patchInfo[thid] = IntervalPad[thidA/RSAA_SAMPLE];
	__syncthreads();

	int patchIndex = (int)( (patchInfo[0]<<(pixelBits+pixelBits))>>(pixelBits+pixelBits) );
	int segI = PatchIndexInverse[patchIndex];
	int segInd[3], segIndd[3];
	segInd[0] = segI/asi.segShift[0];
	segInd[1] = segI%asi.segShift[0]/asi.segShift[1];
	segInd[2] = segI%asi.segShift[1];
	segIndd[0] = segInd[0]+1;
	segIndd[1] = segInd[1]+1;
	segIndd[2] = segInd[2]+1;

	
	extern __shared__ float func[];
	int funcSize = (asi.deg[0]+1)*(asi.deg[1]+1)*(asi.deg[2]+1);

	if(thid<funcSize)
		func[thid] = FuncData[patchIndex*funcSize+thid];
	
	__shared__ uint hit[PIXELTHREADS];
	hit[thid] = 0;
	__syncthreads();

	if(patchInfo[thid]==0&&thid!=0)				{return;}

	int indOff;	
	uint x = (patchInfo[thid]>>(patchBits+pixelBits));
	uint y = (patchInfo[thid]<<pixelBits)>>(patchBits+pixelBits);
	
	float AAoffset[2][2][2]={{{-0.125,-0.375},{0.375,-0.125}},{{-0.375,0.125},{0.125,0.375}}};

	float screen[3] ={2*(x+AAoffset[thid&1][(thid&2)>>1][0]+0.5f)/cRenderParam.renderWidth-1.0f, 1.0f-2*(y+AAoffset[thid&1][(thid&2)>>1][1]+0.5f)/cRenderParam.renderHeight, 0.0f};
	float dir[3];

	float w = screen[0]*cRenderParam.invViewProj[0][3]+screen[1]*cRenderParam.invViewProj[1][3]+screen[2]*cRenderParam.invViewProj[2][3]+cRenderParam.invViewProj[3][3];
	dir[0] = (screen[0]*cRenderParam.invViewProj[0][0]+screen[1]*cRenderParam.invViewProj[1][0]+screen[2]*cRenderParam.invViewProj[2][0]+cRenderParam.invViewProj[3][0])/w-cRenderParam.eye[0]; 
	dir[1] = (screen[0]*cRenderParam.invViewProj[0][1]+screen[1]*cRenderParam.invViewProj[1][1]+screen[2]*cRenderParam.invViewProj[2][1]+cRenderParam.invViewProj[3][1])/w-cRenderParam.eye[1]; 
	dir[2] = (screen[0]*cRenderParam.invViewProj[0][2]+screen[1]*cRenderParam.invViewProj[1][2]+screen[2]*cRenderParam.invViewProj[2][2]+cRenderParam.invViewProj[3][2])/w-cRenderParam.eye[2]; 


	float tmin, tmax;
	RayBoxIntersect(&tmin, &tmax, cRenderParam.eye, dir, asi.box0, asi.boxSize, segInd, segIndd);
	//if(tmax<1||fabsf(tmin)>1.0e5||fabsf(tmax)>1.0e5)	{hit[thid]=0;}
	//if(tmin<1)	tmin=1;

	float fValue[MAX_DEG_INTERPOLATION+1], coeff[MAX_DEG_INTERPOLATION+1];

	//********* code surround by /**/ is chebyshev interpolation ****************//
	//remember: change the interpolation here will also need change the InvMatrix;
	float In[3]={(cRenderParam.eye[0]+tmin*dir[0]-asi.box0[0])/asi.boxSize[0]-segInd[0],
		(cRenderParam.eye[1]+tmin*dir[1]-asi.box0[1])/asi.boxSize[1]-segInd[1], 
		(cRenderParam.eye[2]+tmin*dir[2]-asi.box0[2])/asi.boxSize[2]-segInd[2]};

	float Out[3];
	Out[0] = In[0]+(tmax-tmin)*dir[0]/asi.boxSize[0];
	Out[1] = In[1]+(tmax-tmin)*dir[1]/asi.boxSize[1];
	Out[2] = In[2]+(tmax-tmin)*dir[2]/asi.boxSize[2];

	float xyz[3];
	float cheby;
	for(int i=0; i<=degSum; i++)
	{
		cheby = (1-(__cosf(PI*(2*i+1)/(2*(degSum+1)))))/2;
		xyz[0] = (1-cheby)*In[0]+cheby*Out[0];
		xyz[1] = (1-cheby)*In[1]+cheby*Out[1];
		xyz[2] = (1-cheby)*In[2]+cheby*Out[2];
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], asi.deg[0], asi.deg[1], asi.deg[2], func, fValue+i);
	}
	//********* code surround by /**/ is chebyshev interpolation ****************//

	indOff = 0;
	for(int i=0; i<=degSum; i++)
	{
		coeff[i] = 0.0;
		for(int j=0; j<=degSum; j++)
		{
			coeff[i] += fValue[j]*cInvMatrix[indOff++];
		}
	}
	int indC;
	for(indC=0; indC<degSum; indC++)
	{
		if(coeff[indC]*coeff[indC+1]<0)	break;
	}



	if(indC!=degSum)
//	if(cRenderParam.material.transparent){}
//	else
	{
#if _KI_BEZ_
		float firstRoot = KnotInsertBez(coeff, degSum);
#else
		float firstRoot = KnotInsert(coeff, degSum);
#endif
		if(firstRoot>0)		{hit[thid] = 1;}
	}
	__syncthreads();

	if(thid%RSAA_SAMPLE==0)
	{
		float* pixel = (float*)(surface + y*pitch) + 4*x;
		float* pixelN = (float*)(surfNorm + y*pitch) + 4*x;
		float alpha2 = (hit[thid]+hit[thid+1]+hit[thid+2]+hit[thid+3]+0.1)/RSAA_SAMPLE;
		float alpha1 = 1-alpha2;


		float rgb[3] = {pixel[0], pixel[1], pixel[2]};
		int neighbor[4] = {-4, 4, -pitch/4, pitch/4};
		int nei = (int)(pixelN[3])-1;
		if(nei<0)	return;
		float* pixel2 = pixel+neighbor[nei];
		float rgb2[3] = {pixel2[0], pixel2[1], pixel2[2]};
		pixel[0] = rgb[0]*alpha1+rgb2[0]*alpha2;
		pixel[1] = rgb[1]*alpha1+rgb2[1]*alpha2;
		pixel[2] = rgb[2]*alpha1+rgb2[2]*alpha2;
		pixel[3] = RENDER_AA;
#if SELF_SHADOW
		float* pixelN2 = pixel2+neighbor[nei];
		pixelN[0] = pixelN2[0];
		pixelN[1] = pixelN2[1];
		pixelN[2] = pixelN2[2];
#endif

	}
}
